| title: | |
| author: - Clovis Deletre - Charles Vitry date: output: html_notebook: theme: cerulean number_sections: no toc: yes toc_float: true editor_options: markdown: wrap: 72 |
Nous souhaitons réalisé l’étude d’une série temporelle et faire des prévisions sur celle-ci.
Cette série temporelle est le trafic mensuel d’une Compagnie aérienne de janvier 2011 à août 2019.
Nos prévisions portent sur les 8 mois de l’année 2019
Import de la base, on sélectionne la colonne des valeurs
library(readr)
data <- read_delim("Trafic-voyageurs.csv",
delim = ";", locale = locale(encoding = "ISO-8859-1"))
Rows: 104 Columns: 2
-- Column specification -----------------------------------------------------------------------------------------------------------------------------------------
Delimiter: ";"
chr (1): dates
dbl (1): trafic
i Use `spec()` to retrieve the full column specification for this data.
i Specify the column types or set `show_col_types = FALSE` to quiet this message.
summary(data)
dates trafic
Length:104 Min. :220876
Class :character 1st Qu.:297154
Mode :character Median :355178
Mean :354651
3rd Qu.:407331
Max. :505190
data_value <- data[,2]
Création de la série chronologique :
library(TSstudio)
data_ts <- ts(data_value, start=2011, frequency=12)
plot_1_TimeSeries(data_ts)
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 methods overwritten by 'htmltools':
method from
print.html tools:rstudio
print.shiny.tag tools:rstudio
print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Attachement du package : ‘plotly’
L'objet suivant est masqué depuis ‘package:ggplot2’:
last_plot
L'objet suivant est masqué depuis ‘package:stats’:
filter
L'objet suivant est masqué depuis ‘package:graphics’:
layout
#revoir l affichage car ca prend pas en compte tt 2019
data_ts_train <- window(data_ts, start = c(2011, 1), end = c(2018,12))
data_ts_test <- window(data_ts, start= c(2019,1), end = c(2019,8))
plot(data_ts, xlim=c(2011,2020))
lines(data_ts_test, col=3)
legend("topleft", lty = 1, col=c(1,3), legend=c("Série chronologique Train", "Série chronologique Test"))
-> strong trend -> patern qui se repete, saisonnalité ?
Analyse de la saisonnalité en superposant chaque année (par mois):
-> en supprimant la tendance on voit bien la saisonnalité => saisonnalité régulière
ggseasonplot(data_ts)
data_ts_without_trend = diff(data_ts)
ggseasonplot(data_ts_without_trend)
DECOMPOSITION : additive / Multiplicative Ts = Trend + Seasonal + Random / Ts = Trend * Seasonal * Random
decomposed_data <- decompose(data_ts_train, type="additive")
plot(decomposed_data$trend)
plot(decomposed_data$seasonal)
plot(decomposed_data$random)
boxplot(data_ts ~ cycle(data_ts))
-> on distingue des saisonnalités => faire régression ca n’a pas de sens => modèle de Buys Ballot
-> bonne repartition du bruit -> quelques outliers
checkresiduals(remainder(decomposed_data))
Warning in modeldf.default(object) :
Could not find appropriate degrees of freedom for this model.
On a tendances + saisonnalité
library(forecast)
mean <- meanf(data_ts_train, h=8)
naivem <- naive(data_ts_train, h=8)
driftm <- rwf(data_ts_train, h=8, drif=T)
snaivem <- snaive(data_ts_train, h=8)
plot(mean, plot.conf = F, main="")
lines(naivem$mean, col=2, lty=1)
lines(driftm$mean, col=5, lty=1)
lines(snaivem$mean, col = 4, lty=1)
legend("topleft", lty=1, col=c(1,2,3,4), legend=c("Mean Method", "Naive Method", "Drif Method", "Seasonal Naive"))
#comparaison :
plot(snaivem, plot.conf = F, main="")
lines(data_ts_test, col = 6, lty=1, lwd=3)
plot(driftm, plot.conf = F, main="")
lines(data_ts_test, col = 6, lty=1, lwd=3)
MASE : Mean Absolute Scaled Error : MAPE : Mean Absolute Percentage Error :
res = pred - val MAE = sum(abs(res))/length(val) RSS = sum(res^2) MSE = RSS/length(val) RMSE = sqrt(MSE)
La plus populaire est la MAPE
MAPE(y_pred, y_true)
$MAPE = (1/n) * Σ(|actual – forecast| / |actu0al|) * 10
“a MAPE value of 6% means that the average difference between the forecasted value and the actual value is 6%”
print(summary(mean))
checkresiduals(mean)
accuracy(mean, data_ts_test)
print(summary(naivem))
checkresiduals(naivem)
accuracy(naivem, data_ts_test)
print(summary(driftm))
checkresiduals(driftm)
accuracy(driftm, data_ts_test)
print(summary(snaivem))
checkresiduals(snaivem)
accuracy(snaivem, data_ts_test)
https://mpra.ub.uni-muenchen.de/77718/1/MPRA_paper_77718.pdf page 175
L’approche de BUYS-BALLOT consiste à introduire des variables indicatrices correspondant à chaque saison définit par le cycle d’observation. Pour les données trimestrielles, on intègre 4 variables indicatrices. Et pour les données mensuelles, on intègre 12 variables indicatrices.
Le modèle doit alors être estimé (sans constante) avec ces variables indicatrices.
Préparation des données.
Annees=time(data_ts_train)
ts_DataFrame =data.frame(trafic=data_ts_train,X=as.numeric(Annees))
Création du modèle
Regression <- lm(trafic~X,data = ts_DataFrame)
\(Xt = Zt + St + \mu t\)
La tendance Prédiction sur les données futurs.
tendance=predict(Regression)
Annees=as.numeric(Annees)
AnneeMoisNumeric=seq(max(Annees)+1/12,length=36,by=1/12)
tendance2=predict(Regression, newdata=data.frame(X=AnneeMoisNumeric))
ts_DataFrame$trafic_residual <- residuals(Regression)
Définissons le mois
ts_DataFrame$mois <- round(ts_DataFrame$X - trunc(ts_DataFrame$X),digit=4)
Création du 2nd modèle avec les mois
Regression2 =lm(trafic_residual~0+as.factor(mois),data=ts_DataFrame)
Prédiction de la saisonnalité
prediction2 =predict(Regression2)
Prédiction sur les mois
MoisNumeric= round(AnneeMoisNumeric - trunc(AnneeMoisNumeric
),4)
Prediction3 =predict( Regression2, newdata= data.frame(mois=MoisNumeric))
Calculons une région de confiance avec l’erreur d’ajustement
ResidusRegression2=residuals(Regression2)
hist(ResidusRegression2)
1.96*sqrt(var(ResidusRegression2))
[1] 19226.16
Le bruit est-il auto corrélé ?
plot(acf(ResidusRegression2))
Affichage de la tendance
Buys_ballot_plot_tendance <- plot(data_ts,
main = "Application du modèle de Buys_Ballot",
xlab = "Années",
ylab = "Nombre de Voyageurs")
#droite de tendance
lines(Annees,tendance,col="blue",lwd=2)
#prédiction de la tendance futur
lines(AnneeMoisNumeric,tendance2,col="red")
NA
NA
Affichage du modèle de Buys Ballot
Buys_ballot_plot <- plot(data_ts,
main = "Application du modèle de Buys_Ballot",
xlab = "Années",
ylab = "Nombre de Voyageurs")
#prédiction du modèle de Buys ballot
lines(Annees,tendance+prediction2,col="blue",lwd=2)
#Interval de confiance
polygon(c(AnneeMoisNumeric,rev(AnneeMoisNumeric)),
c(tendance2+Prediction3-1.96*sqrt(var(ResidusRegression2)),
rev(tendance2+Prediction3+1.96*sqrt(var(ResidusRegression2)))),
col="cadetblue1",border=NA)
#Prediction des valeurs
lines(AnneeMoisNumeric,tendance2+Prediction3,col="blue",lwd=2)
lines(data_ts_test,col="black",lwd=3)
Affichage de la prédiction sur les 8 mois de 2020
Buys_ballot_plot <- plot(data_ts_test,
main = "Application du modèle de Buys_Ballot",
xlab = "Années",
ylab = "Nombre de Voyageurs")
#prédiction du modèle de Buys ballot
lines(Annees,tendance+prediction2,col="blue",lwd=2)
#Interval de confiance
polygon(c(AnneeMoisNumeric,rev(AnneeMoisNumeric)),
c(tendance2+Prediction3-1.96*sqrt(var(ResidusRegression2)),
rev(tendance2+Prediction3+1.96*sqrt(var(ResidusRegression2)))),
col="cadetblue1",border=NA)
#Prediction des valeurs
lines(AnneeMoisNumeric,tendance2+Prediction3,col="blue",lwd=2)
lines(data_ts_test,col="black",lwd=3)
Nous avons réussi à ajuster une droite de régression. on remarque que la prédiction semble bien correspondre à la réalité si on fait abstraction du dernier mois où le nombre de voyageurs a bien plus chuté que la prédiction du modèle de Buys-Balot.
Comparons avec un ajustement local réalisé par lissage moyennes mobiles.
Mettre belle formule en latex ici
fcst_se <- ses(data_ts_train, h = 8)
print(summary(fcst_se))
Forecast method: Simple exponential smoothing
Model Information:
Simple exponential smoothing
Call:
ses(y = data_ts_train, h = 8)
Smoothing parameters:
alpha = 0.2559
Initial states:
l = 258126.0245
sigma: 31480.96
AIC AICc BIC
2430.727 2430.988 2438.420
Error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 7057.127 31151.31 25752.89 1.326234 7.684321 1.002462 0.0711143
Forecasts:
checkresiduals(fcst_se)
Ljung-Box test
data: Residuals from Simple exponential smoothing
Q* = 144.66, df = 17, p-value < 2.2e-16
Model df: 2. Total lags used: 19
plot(fcst_se)
lines(data_ts_test, col="red")
df_se = as.data.frame(fcst_se)
predict_value_se <- df_se$`Point Forecast`
MAPE(predict_value_se, data_ts_test)*100
[1] 7.658334
Fit Exponential Smoothing model -> trouve le meilleur lissage expo
fit_ets <- ets(data_ts_train)
print(summary(fit_ets))
ETS(A,A,A)
Call:
ets(y = data_ts_train)
Smoothing parameters:
alpha = 0.1568
beta = 1e-04
gamma = 1e-04
Initial states:
l = 248267.1099
b = 2163.3982
s = -17928.3 -29535.73 9295.935 11005.81 -57117.85 -7708.17
38272.64 14592.34 16899.53 34763.15 -7344.204 -5195.15
sigma: 11014.45
AIC AICc BIC
2241.611 2249.458 2285.205
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -458.6799 10054.77 7831.554 -0.2623253 2.371375 0.3048527 0.09626331
checkresiduals(fit_ets)
Ljung-Box test
data: Residuals from ETS(A,A,A)
Q* = 7.1794, df = 3, p-value = 0.06639
Model df: 16. Total lags used: 19
fcst_ets <- forecast(fit_ets, h=8)
plot(fcst_ets)
lines(data_ts_test, col="red")
df_ets = as.data.frame(fcst_ets)
predict_value_ets = df_ets$`Point Forecast`
MAPE(predict_value_ets, data_ts_test)*100
[1] 3.005848
# retourne les meilleurs paramètres
# d=1 enleve la tendance
# D=1 enleve la saisonnalité
# => avoir des données stationnaires
# trace : voir les résultats
fit_arima <- auto.arima(data_ts_train, d=1, D=1, stepwise = FALSE, approximation = FALSE, trace=TRUE)
ARIMA(0,1,0)(0,1,0)[12] : 1846.398
ARIMA(0,1,0)(0,1,1)[12] : 1833.134
ARIMA(0,1,0)(0,1,2)[12] : 1835.211
ARIMA(0,1,0)(1,1,0)[12] : 1833.056
ARIMA(0,1,0)(1,1,1)[12] : 1835.09
ARIMA(0,1,0)(1,1,2)[12] : Inf
ARIMA(0,1,0)(2,1,0)[12] : 1835.207
ARIMA(0,1,0)(2,1,1)[12] : 1837.012
ARIMA(0,1,0)(2,1,2)[12] : 1836.461
ARIMA(0,1,1)(0,1,0)[12] : 1814.951
ARIMA(0,1,1)(0,1,1)[12] : 1801.155
ARIMA(0,1,1)(0,1,2)[12] : 1803.362
ARIMA(0,1,1)(1,1,0)[12] : 1803.592
ARIMA(0,1,1)(1,1,1)[12] : 1803.361
ARIMA(0,1,1)(1,1,2)[12] : Inf
ARIMA(0,1,1)(2,1,0)[12] : 1805.004
ARIMA(0,1,1)(2,1,1)[12] : 1805.397
ARIMA(0,1,1)(2,1,2)[12] : Inf
ARIMA(0,1,2)(0,1,0)[12] : 1816.915
ARIMA(0,1,2)(0,1,1)[12] : 1803.033
ARIMA(0,1,2)(0,1,2)[12] : 1805.296
ARIMA(0,1,2)(1,1,0)[12] : 1805.702
ARIMA(0,1,2)(1,1,1)[12] : 1805.295
ARIMA(0,1,2)(1,1,2)[12] : Inf
ARIMA(0,1,2)(2,1,0)[12] : 1807.026
ARIMA(0,1,2)(2,1,1)[12] : 1807.441
ARIMA(0,1,3)(0,1,0)[12] : 1817.787
ARIMA(0,1,3)(0,1,1)[12] : Inf
ARIMA(0,1,3)(0,1,2)[12] : Inf
ARIMA(0,1,3)(1,1,0)[12] : Inf
ARIMA(0,1,3)(1,1,1)[12] : Inf
ARIMA(0,1,3)(2,1,0)[12] : Inf
ARIMA(0,1,4)(0,1,0)[12] : 1820.052
ARIMA(0,1,4)(0,1,1)[12] : Inf
ARIMA(0,1,4)(1,1,0)[12] : Inf
ARIMA(0,1,5)(0,1,0)[12] : Inf
ARIMA(1,1,0)(0,1,0)[12] : 1825.579
ARIMA(1,1,0)(0,1,1)[12] : 1812.512
ARIMA(1,1,0)(0,1,2)[12] : 1814.657
ARIMA(1,1,0)(1,1,0)[12] : 1813.2
ARIMA(1,1,0)(1,1,1)[12] : 1814.614
ARIMA(1,1,0)(1,1,2)[12] : 1816.192
ARIMA(1,1,0)(2,1,0)[12] : 1815.227
ARIMA(1,1,0)(2,1,1)[12] : Inf
ARIMA(1,1,0)(2,1,2)[12] : 1817.796
ARIMA(1,1,1)(0,1,0)[12] : 1816.841
ARIMA(1,1,1)(0,1,1)[12] : 1802.853
ARIMA(1,1,1)(0,1,2)[12] : 1805.117
ARIMA(1,1,1)(1,1,0)[12] : 1805.653
ARIMA(1,1,1)(1,1,1)[12] : Inf
ARIMA(1,1,1)(1,1,2)[12] : Inf
ARIMA(1,1,1)(2,1,0)[12] : Inf
ARIMA(1,1,1)(2,1,1)[12] : Inf
ARIMA(1,1,2)(0,1,0)[12] : 1819.234
ARIMA(1,1,2)(0,1,1)[12] : 1805.22
ARIMA(1,1,2)(0,1,2)[12] : 1807.539
ARIMA(1,1,2)(1,1,0)[12] : 1807.381
ARIMA(1,1,2)(1,1,1)[12] : 1807.538
ARIMA(1,1,2)(2,1,0)[12] : 1808.925
ARIMA(1,1,3)(0,1,0)[12] : 1820.05
ARIMA(1,1,3)(0,1,1)[12] : 1806.055
ARIMA(1,1,3)(1,1,0)[12] : 1808.732
ARIMA(1,1,4)(0,1,0)[12] : Inf
ARIMA(2,1,0)(0,1,0)[12] : 1824.435
ARIMA(2,1,0)(0,1,1)[12] : 1811.07
ARIMA(2,1,0)(0,1,2)[12] : 1813.287
ARIMA(2,1,0)(1,1,0)[12] : 1811.619
ARIMA(2,1,0)(1,1,1)[12] : 1813.247
ARIMA(2,1,0)(1,1,2)[12] : Inf
ARIMA(2,1,0)(2,1,0)[12] : 1813.821
ARIMA(2,1,0)(2,1,1)[12] : 1815.872
ARIMA(2,1,1)(0,1,0)[12] : Inf
ARIMA(2,1,1)(0,1,1)[12] : Inf
ARIMA(2,1,1)(0,1,2)[12] : Inf
ARIMA(2,1,1)(1,1,0)[12] : Inf
ARIMA(2,1,1)(1,1,1)[12] : Inf
ARIMA(2,1,1)(2,1,0)[12] : Inf
ARIMA(2,1,2)(0,1,0)[12] : Inf
ARIMA(2,1,2)(0,1,1)[12] : Inf
ARIMA(2,1,2)(1,1,0)[12] : Inf
ARIMA(2,1,3)(0,1,0)[12] : Inf
ARIMA(3,1,0)(0,1,0)[12] : 1823.646
ARIMA(3,1,0)(0,1,1)[12] : 1808.49
ARIMA(3,1,0)(0,1,2)[12] : 1810.542
ARIMA(3,1,0)(1,1,0)[12] : 1808.594
ARIMA(3,1,0)(1,1,1)[12] : 1810.321
ARIMA(3,1,0)(2,1,0)[12] : 1810.708
ARIMA(3,1,1)(0,1,0)[12] : Inf
ARIMA(3,1,1)(0,1,1)[12] : Inf
ARIMA(3,1,1)(1,1,0)[12] : Inf
ARIMA(3,1,2)(0,1,0)[12] : Inf
ARIMA(4,1,0)(0,1,0)[12] : 1823.996
ARIMA(4,1,0)(0,1,1)[12] : 1810.199
ARIMA(4,1,0)(1,1,0)[12] : 1810.845
ARIMA(4,1,1)(0,1,0)[12] : Inf
ARIMA(5,1,0)(0,1,0)[12] : 1825.055
Best model: ARIMA(0,1,1)(0,1,1)[12]
print(summary(fit_arima))
Series: data_ts_train
ARIMA(0,1,1)(0,1,1)[12]
Coefficients:
ma1 sma1
-0.7675 -0.5465
s.e. 0.0977 0.1295
sigma^2 = 138827719: log likelihood = -897.43
AIC=1800.85 AICc=1801.16 BIC=1808.11
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 805.2378 10822.93 7747.841 0.2158443 2.213481 0.3015941 0.03453594
checkresiduals(fit_arima)
Ljung-Box test
data: Residuals from ARIMA(0,1,1)(0,1,1)[12]
Q* = 10.898, df = 17, p-value = 0.8618
Model df: 2. Total lags used: 19
fcst_arima <- forecast(fit_arima, h=8)
plot(fcst_arima)
lines(data_ts_test, col='red')
df_arima = as.data.frame(fcst_arima)
predict_value_arima = df_arima$`Point Forecast`
MAPE(predict_value_arima, data_ts_test)*100
[1] 2.814135